Critical comparison of electrode models in density functional theory based quantum 

transport calculations 
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We study the performance of two different electrode models in quantum transport calculations 
based on density functional theory: Parametrized Bethe lattices and quasi-one dimensional wires or 
nanowires. A detailed account of implementation details in both cases is given. From the systematic 
study of nanocontacts made of representative metallic elements, we can conclude that parametrized 
electrode models represent an excellent compromise between computational cost and electronic 
structure definition as long as the aim is to compare with experiments where the precise atomic 
structure of the electrodes is not relevant or defined with precision. The results obtained using 
parametrized Bethe lattices are essentially similar to the ones obtained with quasi one dimensional 
electrodes for large enough sections of these, adding a natural smearing to the transmission curves 
that mimics the true nature of polycrystalline electrodes. The latter are more demanding from the 
computational point of view, but present the advantage of expanding the range of applicability of 
transport calculations to situations where the electrodes have a well-defined atomic structure, as is 
case for carbon nanotubes, graphene nanoribbons or semiconducting nanowires. All the analysis is 
done with the help of codes developed by the authors which can be found in the quantum transport 
toolbox Alacant and are publicly available. 



I. INTRODUCTION 

One of the most active research fields in Nanoscience 
is the one focusing on understanding and controlling the 
charge transport between bulk electrodes when these are 
connected by an atomic- or a molecular-size region and a 
bias voltage is applied between them[l[, in other words, 
understanding and controlling the formation and con- 
comittant resistance of nanoscopic bridges. More than 
10 years of experimental along with theoretical work is 
finally taking us to an unprecedented level of control 
and comprehension of these systems. On the theoreti- 
cal side, we have witnessed the marriage of theoretical 
quantum transport basics and density functional theory, 
giving birth to one of the most active and fructiferous 
fields in theoretical nanoscience (2|-[68| . 

For nanoscopic conductors every atom counts and the 
transport properties are strongly dependent on the de- 
tailed atomic arrangement. Thus, in order to make 
theoretical predictions that can be compared with ex- 
perimental results, it is important, to have a reliable 
description of, first, the atomic structure of the con- 
ductor and, second, the accompanying electronic struc- 
ture. This can be achieved most conveniently with the 
aid of ab initio electronic structure methods based on 
atomic orbitals such as, e.g., Gaussian^, Crystal [70l| 
or Siesta[71[. These codes implement density func- 
tional theory (DFT) to obtain an effective mean-field de- 
scription of the electronic structure of, in our case, the 
nanoscopic bridge. This is typically done through the 
effective Kohn-Sham one-body Hamiltonian that takes 
into account the electron-electron interactions at a static 
mean-field level. 

A central challenge in the theoretical description of 
these systems is that the electronic structure of the 



atomic- or molecular-size conductor is altered by the cou- 
pling to the bulk electrodes. Thus, in calculating the elec- 
tronic structure the coupling to the (semi-infinite) elec- 
trodes has to be taken into account. This poses the dif- 
ficult problem of dealing with an infinite system without 
translation invariance. In addition, one should strictly 
carry out the electronic and atomic structure calcula- 
tion out of equilibirum, as imposed by the applied volt- 
age. This is usually done by making use of the one-body 
Green's functions (GFs) and the so-called partitioning 
technique, as will be explained in the following sections. 

Clearly, while the detailed atomic and electronic struc- 
ture of the nanoscopic bridge plays a crucial role, farther 
away from the bridge these become less important. Be- 
sides, the exact atomic structure of the bulk electrodes, 
as encountered in real experiments, is not known and 
cannot be controlled with precision beyond a few contact 
atoms [7^|. This lack of control lies behind the statisti- 
cal deviations observed in measurable quantitites such as 
the conductance. This brings us to the important ques- 
tion of how to model the bulk electrodes, mantaining a 
compromise between realism and computational effort. 
Several possibilities of how to model the electrodes have 
been presented in the literature, with every model having 
advantages and disadvantages [1, [E fl • Here we explore 
the use of two types of electrode models: (i) parametrized 
tight-binding (TB) Bethe lattices (fj and (ii) perfect 
nanowires of finite section, stressing their weaknesses and 
strengths as models to represent the reality. Another re- 
lated question which we partially address in this paper is 
to what extent the particular shape or atomic arrange- 
ment of the electrodes near the bridge introduces varia- 
tions in the conductance and how these depend on the 
chemical nature of the atoms. The results presented be- 
low are all obtained with codes developed by the authors 
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FIG. 1: Sketch of the transport problem. The system is 
divided into three parts: Left electrode (L), device (D), and 
right electrode (R). 



over the years which can be found in the publicly avail- 
able quantum transport toolbox ALACANT[73l]. 



II. NON-EQUILIBRIUM GREEN'S FUNCTIONS 
AND LANDAUER FORMALISMS 

In the following, and for completeness' sake, we give 
a summary of the central aspects to the non-equilibrium 
Green's function (NEGF) and Landauer formalisms for- 
mulated for a non-orthogonal localized atomic basis set. 
Although most of the details can be found in the early 
literature \Ml MM, M, El [H], here we discuss 
in depth some of those that are not usually addressed 
and become essential for a correct implementation of the 
above mentioned formalisms. 

We divide the system into three parts, as shown in Fig. 
[TJ The semi-infinite left (L) and right (R) electrodes or, 
hereon, leads, and the intermediate region between the 
two leads hereon called device (D) which contains the 
a central, narrow region where most of the scattering 
takes place (e.g., a nanoscopic constriction of the same 
material as the leads or a trapped molecule) . We assume 
that the leads are only coupled to the scattering region 
but not to each other. In a localized atomic basis set the 
Hamiltonian H of the system is given by: 

/ H L H LD \ 
H = H DL H D H DR (1) 
\ H RD H r J 

Since atomic basis sets are usually non-orthogonal we 
also have to take into account the overlap between the 
atomic orbitals given by the matrix S: 

/ s L S LD \ 
S = S DL S D S DR (2) 
\ S RD S R J 

In order to deal with the problem of an infinite system 
without translation invariance, it is convenient to make 
use of the one-body Green's functions as explained, e.g., 
in the book by Economou(74|. The one-body Green's 
function (GF) is defined as the resolvent operator of the 



one-body Schrodingcr equation: 

(z - H)G(z) = I. (3) 

where z is, in general, a complex number and I is the 
identity. 

If z does not coincide with an eigenvalue tk of the 
Hamiltonian H the GF operator has the following simple 
solution: G(z) = (z — H)^ 1 for z^ft. 

Obviously, for z — the GF operator has a pole and 
is thus not well defined. In this case one can define two 
GFs by a limiting procedure: The regarded GF is defined 
as & + '(E) := lim^o G(E + irj), and the advanced GF is 
defined as &~\E) := lim^o G(E—irj) where E is a real 
number (the energy). The retarded (advanced) GF can 
be analytically continued into the upper (lower) complex 
plane. Moreover, away from the poles of G(z), i.e., for 
z =^ eft both definitions coincide with G(z): & + \z) = 
G(-)(z) = G{z). 

Because of the non-orthogonality of the basis set it is 
convenient to define the following Green's function ma- 
trix 

(zS - H)G(z) = 1. (4) 

Note that this Green's function matrix is not the stan- 
dard one which is simply defined by the matrix elements 
of the GF operator Gr(,z) j). However, the latter GF 
matrix is inconvenient to handle in the case of non- 
orthogonal basis sets (see Appendix [S] for a complete 
discussion). 

Using the technique explained in App. [B] it can be 
shown that the GF of the device region is given by the 
following matrix: 

G D (^) = (zS D -H D -S L (z)-S R (z))- 1 (5) 

where Sl and are the so-called lead self-energies 
which describe the coupling of the device to the semi- 
infinite L and R leads. These self-energies can be calcu- 
lated from the Green's functions of the (isolated) leads, 
g a (z) = (zS a -H a )-h 

E a («) = (zS DQ - H DQ ) Sa (z) (zSl a - H+ a ) (6) 

where the index a denotes the electrode L or R, and we 
have exploited the hermiticity of the Hamiltonian and 
the overlap matrix, i.e., = H q d and S^ Da = S a u- 

All quantities of interest such as the density of states 
(DOS), charge density, current /, and zero-bias as well 
as differential conductance dl/dV can be calculated from 
the GF matrix of the device region Gd and the lead self- 
energies Sl and Sr. For instance, in the case of an 
effective Kohn-Sham one-body Hamiltonian, as consid- 
ered here, the current through the nanoscopic conductor 
is given by the famous Landauer formula [75l|: 

I = ^JdE [f(E - p L ) - f(E - Mr)] T(E), (7) 
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where / represents the Fermi distribution function, /i Q 
the left and right chemical potentials, and the transmis- 
sion function, T(E), can be calculated from the retarded 
and advanced GFs by the Caroli expression (76j: 



T(E) = Tr ^Tl(E)G^\e)T-r(E)G) [ ^ > (£7) J . (8) 

Here IY and Tr are the so-called coupling matrices which 
are defined as 



(+)/ 



T h (E) := i (s[ + )(£)-S[->(£) 



.(+) 



(E) 



-(-) 



(E) 



(9) 
(10) 



Note that, since the self-energy matrices are usually sym- 
metric, the coupling matrices are just (twice) the imagi- 
nary parts of the self-energies. 

At zero temperature the zero-bias conductance is now 
just given by the transmission function at the Fermi level 
fi (i.e. the electrochemical potential at zero tempera- 
ture): 



G 



2e 2 



xT{fi) 



(11) 



Hence, the transmission function is the central quan- 
tity for calculating the electronic transport properties of 
nanoscopic conductors. It is worth noting at this point 
that there is a controversy on the use of (Kohn-Sham) 
DFT to calculate the transmission function. In addition 
to the obvious fact that there is no mathematical sup- 
port to the use of DFT out of equilibrium, it has been 
recently argued that a DFT description of the device re- 
gion can never yield the right value of the zero-bias con- 
ductance, not even using the exact exchange-correlation 
potential in case this was known. In fact, the correc- 
tions to the DFT zero-bias transmission calculated using 
standard functionals can be important in high resistance 
cases. We refer the interested reader to Ref. [73] for a 
full discussion of these issues which are beyond the scope 
of this work. 



III. ELECTRODE MODELS 

In the ALACANT toolbox two different codes can be 
found, differing in the way the bulk electrodes are im- 
plemented. In the first the electrodes can be described 
by a parametrized TB Bethc lattice (BL) model with the 
coordination and parameters appropriate for the chosen 
electrode material. In the second the electrodes are ap- 
proximated by finite section wires described with a Kohn- 
Sham Hamiltonian, usually computed at the same level 
as that of the scattering region or device. 



A. Bethe Lattice electrodes 

A Bethc lattice, sometimes also called Caley tree, is 
generated by connecting a site with TV nearest-neighbors 




FIG. 2: Left: Finite section of Bethe lattice (BL) with coor- 
dination 6. All atoms of the BL have the same coordination as 
in the corresponding crystalline structure giving rise to short 
range order. But there is no long range order in the BL due to 
the absence of closed loops. Right: 2D cartoon of a nanocon- 
tact (big grey circles) with the first atoms of the BL (small 
white circles) attached to the outer planes of the nanocontact. 



in directions that can be those of a particular crystalline 
lattice. The new TV sites are each one of them connected 
to TV — 1 different sites and so on and so forth. The gen- 
crated lattice has the local topology of an actual lattice 
(number of neighbors and crystal directions) but has no 
rings, and thus does not describe the long range order 
characteristic of real crystals. The left hand side of Fig. 
[2] shows the first three layers of a BL with coordination 
6. 

The advantage of choosing a BL over other models 
resides on the one hand in the lack of long-range order 
which mimics the poly-crystallinity of real electrodes. On 
the other hand the BL captures the short-range order 
since the local coordination of an atom is that of an atom 
in the bulk crystal of the corresponding material. 

The right hand side of Figure [5] illustrates schemati- 
cally how the device (represented here by a single-element 
nanocontact) is connected to the BL electrodes: For a 
given chosen atom, typically in the outer planes of the 
device, a branch of the BL is added in the direction r, 
of any missing bulk atom (including those missing in the 
same plane). The directions in which tree branches are 
added are indicated by white small circles which repre- 
sent the first atoms of the branch in that direction. This 
corresponds to adding a BL self-energy S Ti to the atom 
in that direction (see App. [C] for a derivation of this 
formula) : 



H Ti (E) = H Ts [EI - H (S T (£) - ^(E))}- 1 

(12) 

where St is the sum over all self-energies in all directions 
and the bar in f, indicates the opposite direction of Tj, i.e. 
f, = — Tj. The electrode self-energies X L and are thus 
obtained by summing up the directional BL self-energies 
S Ti in all directions Tj missing on that particular atom 
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FIG. 3: Comparison of Bethe Lattice DOS resulting from a parametrization obtained from an LDA Hamiltonian (green dashed 
lines) and from a Papaconstantopoulus parametrization (blue dotted lines) with the DOS for a real FCC lattice calculated in 
LDA (red continuous lines) for the three electrode materials Cu (a), Al (b) and Ni (c). 



S OL/RiTi for all the atoms connected to that electrode: 



£ L/R (£) 



E 

all atoms a. 



E 



j^ /R all missing 
connected to L/R directions Tj 



(13) 

Assuming that the most important structural details of 
the electrode are already included in the central cluster, 
the BLs should have no other relevance than that of in- 
troducing a generic bulk electrode for a given metal. 

In order to compare the results of the BL model with 
the actual electronic structure of the corresponding real 
crystal lattice, we calculate the bulk DOS of the BL from 
the imaginary part of the local GF: 



Po(E) = --lmTr[G {E)}, 

7T 



(14) 



where the local Green's function Go of the Bethe lattice 
is given by 



G (£) = CE-Ho-Er(£))" 



(15) 



In Fig. [3] we compare the bulk DOS of BL models us- 
ing different parametrizations with electronic structure 
calculations in the local density approximation (LDA) 
for the three different electrode materials considered here 
(Cu, Al and Ni). The BLs have coordination 12, corre- 
sponding to the FCC crystalline structure of the bulk 
materials. On the one hand we have taken the TB pa- 
rameters directly from the nearest-neighbor hoppings and 
on-site energies of the LDA Kohn-Sham Hamiltonian of 
the FCC crystal (ignoring the overlap), where the cal- 
culations have been carried out with the help of the 
CRYSTAL code. On the other hand we have taken the 
TB parameters established by Papaconstantopoulus and 
coworkers [Z^. As can also be seen from Fig. [3J the BL 
construction results in all cases in a smooth DOS which 
reproduces the basic features of the one corresponding 
to a mono-crystalline solid. As can be seen in Fig. [3l 
depending on the type of material, the use of one set of 
parameters or another can result in a more accurate de- 
scription. However, the choice of parameters should not 



be determinant in the final conductance results as long 
as the device is sufficiently large. 

As usual we have assumed an orthonormal basis set 
for the Bethe lattice. On the other hand, the basis set 
of the device region is typically non-orthogonal. Hence, 
the question arises of how to match the two different 
levels of modeling. A straightforward approach is to or- 
thogonalize the device basis set, for example with the 
Lowdin orthogonalization scheme through the transfor- 
mation Hd = Sd~ 1//2 HdSd _1 ^ 2 - Equivalently, one can 



J L/R 



■-D - J D 

de-orthogonalize the self-energies Sl and Sr.: S[ 

o 1/2 V q 1/2 

Alternatively, one can obtain the TB parameters in a 
non-orthogonal basis set, and take into account the over- 
lap between orbitals on neighbouring atoms in the cal- 
culation of the BL self-energies. In this case the Dyson 
equation for the calculation of the BL self-energy is triv- 
ially modified as follows: 

H Ti (E) = (H Ti -ES Ti ) (16) 
x [ES Q H (H T (E) - -ErAEW 1 (H£ - ESl t ) 

where S Ti is the overlap matrix between orbitals on neigh- 
boring atoms in the direction Tj. 

In case of a non-orthogonal basis set one has to take 
into account the non-diagonal part of the GF between 
different sites of the BL when computing the BL DOS: 



Po(E) 



1. 



-ImTr 



G (E)S + J2G , n {E)S Ti 



(17) 



This is most easily done by extending the unit cell Hamil- 
tonian of the BL with all nearest neighbour sites, com- 
puting the GF Gxo(E) of the extended unit cell (XO), 
and then taking the partial trace for the central site of 
the matrix product Gxo(E)Sxo- 

Po (E) = --ImTr [G X o(£)Sxo]. (18) 

In Fig. 2]we compare the bulk DOS of BL models with 
different parametrizations (taking into account the over- 
lap between orbitals on neighbouring atoms) with LDA 
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FIG. 4: Bethe Lattice DOS with different parametriztions 
taking into account overlap between atomic orbitals com- 
pared to DOS of real FCC lattice calculated with LDA (red 
continuous lines) for Al. The BL parameters have been 
taken either directly from the LDA Kohn-Sham Hamiltonian 
(green dashed lines) or from the Papaconstantopoulus TB 
parametrization with overlap (blue dotted lines). 




FIG. 5: Sketch of transport problem for the case of one- 
dimensional nanowires as electrodes. The system is divided 
into 3 parts: Left electrode (L), device (D) and right electrode 
(R). 



electronic structure calculations for the case of Al. As 
before we have taken the TB parameters either directly 
from the nearest-neighbor hoppings and on-site energies 
of the LDA Kohn-Sham Hamiltonian of the FCC crystal 
(this time taking into account the overlap) or we have 
taken the TB parameters established by Papaconstan- 
topoulus and coworkers (this time with overlap). Now 
the BL DOS does not resemble the DOS of a real crys- 
tal lattice anymore: The band width now becomes semi- 
infinite extending infinitely to positive energies. This ar- 
tifact can only be healed by scaling down the overlap 
considerably. We therefore conclude that the introduc- 
tion of non-orthogonal basis sets in the description of the 
BL does not seem to be very useful for mimicking the 
DOS of real materials, being preferable the self-energy 
deorthogonalization procedure described above. 



B. Nanowire electrodes 

The second type of model for the leads consists of 
finite-section wires where the electronic structure is de- 
scribed at the same computational level as that of the 



device. As indicated in Fig. [5l we subdivide the one- 
dimensional leads into unit cells which must be cho- 
sen sufficiently large so that the coupling between non- 
neighboring unit cells can be neglected. In general a unit 
cell consists of several primitive unit cells of the crys- 
tal. The Hamiltonian matrix of the left lead Hl can be 
subdivided into sub-matrices in the following manner: 



H, 



\ 



H{ H Hr 



V o 



HI H Hr 
Hi Ho/ 



(19) 



Analogously the Hamiltonian of the right lead is given 
by the following matrix: 



H 



/H Hr 
H\ H Hx 



R — 



H| H Hx 



V o 



/ 



(20) 



In a similar way, the overlap inside the leads is given by 
the matrices 



St = 



Sq Si 



Si So Si 



V o 



(21) 



SI So/ 



and 



/So Si a \ 

S{ So Si 
S{ So Si 

V o •• ••• •./ 



(22) 



Furthermore, the unit cell of each lead that is immedi- 
ately connected to the scattering region (unit cell "Z" for 
the left and unit cell "r" for the right lead) is included 
into the device part of the system. So the Hamiltonian 
of the device region reads 



Hr 




H; : s 0/, r 

Hs He, r 
H r s H r 



(23) 



and the overlap matrix is given by: 



s« 


S/,s 


0l, r 


Ss,i 


Ss 


Ss.r 


r ,l 


Sr.S 


s r 



(24) 



Since only the Z- and r-parts of the device region are 
immediately connected to the two semi-infinite nanowire 
electrodes the self-energy matrices Sl and Sr. that de- 
scribe the coupling of the device region to the electrodes 
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L and R are given by matrices that arc different from 
zero only in the I- and r-parts, respectively: 

/ Ej(z) Ol.s Ol,r \ 

X L {z) = S ,l Os s , r (25) 

\ Or.l Or,S Or J 

and 

/ ; 1<S l>r \ 
S R (z) = Os,i Os s , r (26) 

\ r , ; r S V r (z) J 

As shown in App. [D] the non-zero submatriccs £/ and 
S r can be calculated from the Hamiltonian and over- 
lap sub-matrices of the two leads by the following Dyson 
equations: 

£,(*) = ( H I -zSlHzSo-Ho-E^))" 1 

x(H x -zSi) (27) 

H r (z) = (Hx-zSO^So-Hq-S^z))- 1 

x(H\-zS{). (28) 

IV. SELF-CONSISTENT ELECTRONIC 
STRUCTURE OF THE DEVICE 

In the context of standard DFT electronic struc- 
ture calculations of finite or periodic systems the self- 
consistent mean-field electronic potential and the density 
matrix (or Kohn-Sham wave functions) are determined 
by the sole input of the atomic structure of the cluster or 
cell, through the chosen exchange-correlation functional. 
In the context of quantum transport, where the systems 
are infinite, but present no translational invariance, the 
"boundary conditions" imposed by the electrodes play an 
additional and important role. The electronic structure 
of the device region also depends on the model chosen to 
represent the electrodes and the details of how to carry 
out the self-consistency may vary, particularly when out 
of equilibrium. We discuss now two different alterna- 
tives: The embedded cluster approach, associated to the 
use of BLs (see Sec. IIII Ap . and the supercell approach, 
where the electrodes are described by nanowires (see Sec. 

lmB|) . 



A. Embedded cluster approach 

In the embedded cluster approach the electronic struc- 
ture of the infinite system is calculated self-consistently 
only within a finite-size region -the scattering or de- 
vice region containing the nanoconstriction or molecule- 
while the electronic structure of the rest of the system 
(i.e., the two bulk electrodes) is fixed from the very be- 
ginning to that of a simplified parametrized BL model 



(see Sec. IIII A|) . The basic premise here is to set up a de- 
vice region sufficiently large. In other words, a sufficiently 
wide section of the bulk electrodes must be included in 
the device region so that the interface resistance between 
the BL and the device, the former being described at a 
lowest level of approximation than the latter, does not 
contribute significantly to the overall resistance which is 
essentially determined by the intrinsic resistance of the 
device. 

In equilibrium the two leads or electrodes must have 
the same electrochemical potential. If the leads are made 
of different materials with different work functions, an 
overall charge transfer must occur somewhere. This gives 
rise to an electric field, shifting the band structures of the 
two leads relative to each other, and subsequently align- 
ing the electrochemical potentials of the two leads. The 
net effect of the charge transfer on the electronic struc- 
ture of the bulk electrode material outside the device 
can be taken into account by simply shifting the electro- 
chemical potentials (and of course the band structure) of 
the two materials to a common electrochemical poten- 
tial. Leads of the same material but presenting different 
crystallographic orientations might also present different 
work functions, but the BL model cannot account for 
this difference. Notice that we have refrained from be- 
ing too specific about where the charge transfer takes 
place. A localized charge transfer in the device region, 
typically a one- or quasi-one-dimensional system, cannot 
be entirely responsible for the electrochemical alignment 
far away from the device for obvious electrostatic rea- 
sons. One must be cautious with the extent of the re- 
gion necessary for this transfer to take place. Only for 
infinite two-dimensional interfaces between different ma- 
terials the charge accumulates strictly at the interface. 
In the opposite limit of one-dimensional systems in point 
contact, the charge transfer must extend logarithmically 
into the bulk[79j]. Regardless of where the charge transfer 
takes place, anyway, thcrmodynamical equilibrium must 
be reached. 

Taking a common electrochemical potential /i, the rigid 
electrostatic shifts are Al = fJb — an d Ar = n — fj,^ 
where ^ and /x R are the electrochemical potentials (or 
work functions) of the materials of the left and right elec- 
trode, respectively. The lead Hamiltonians corrected by 
the electrostatic shift are thus given by Hl + AlSl and 
Hr + ArSr . Now the Kohn-Sham Hamiltonian of the 
entire system (leads+device) is given by 

H KS - = (29) 

/ Hl-z^Sl H ld -m l Sld \ 

Hdl — Ml^dl Hd — mSd Hdr — Mr^dr 
\ Hrd-MrSrd Hr-/XrSr / 

where we have included the electrochemical potential fi of 
the entire system. The corresponding overlap matrix that 
takes into account the non-orthogonality of the atomic 
orbitals has of course the same form as in eq. [2] 

The Kohn-Sham Hamiltonian of the entire system de- 
pends only on the electron density nu(r) of the device 
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FIG. 6: Diagram illustrating the self-consistent procedure for 
KS based NEGF formalism as explained in the text. 

region since the electronic structure of the rest of the 
system is kept fixed (apart from the electrostatic shifts 
Al and Ar which depend on the electrochemical poten- 
tial n): Hks = Hks["d]- The electron density nn(r) 
can be obtained from the density matrix of the device 
region: [8(| 

nu(f f )= ]T <p a (r)P aP ^(r) (30) 

a,/3£D 

which in turn is found by integrating (most conveniently 
done by analytic continuation to the complex plane) the 
device part of the Green's function: 

PdM = --Im / dEG$\E) (31) 

where the device Green's function is given by 

G [ +\E) = (ES D - H D - £ L (£) - ^r(E))- 1 (32) 

Now in order to obtain the electrochemical potential 
of the entire system, one has to impose charge neutrality 
in the entire system. Since the metallic leads outside the 
device region are charge neutral themselves, it suffices to 
impose charge neutrality within the device region: 

JV d (aO = Tr[P D (/i)S D ] 

= Im / d£Tr[G D +) (#) S D ] (33) 

T J-oo 

Since Gd via Hd is a functional of the electron density 
riD(r), the electronic structure of the device region can 
now be determined sclf-consistcntly, already having taken 
into consideration the effect of the leads through the self- 
energies in Eq. I32I For a practical implementation of this 
procedure we have created an interface to the quantum 
chemistry code Gaussian03, taking thus advantage of 
the various DFT implementations that can be found in 



such a code. Figure [6] shows a schematic picture of the 
self-consistent calculation of the electronic structure of 
the device region as described above. 

Out of equilibrium, i.e., for a finite bias eV = fJ-L — fJ-R 
one has to use the NEGF technique. In this case there 
is an additional contribution to the density matrix which 
can be calculated by integrating the so-called lesser GF 
G < within the bias window: 

PdVl, WO = PdVr) - f dEG<(E), (34) 

where we have assumed a positive difference between L 
and R electrochemical potentials, hl — Hr, > 0. As for the 
equilibrium case, charge neutrality must be imposed, e.g., 
by setting /ir, to the appropriate value. The lesser GF 
can be easily calculated from the retarded and advanced 
GFs: 

G<(E) = iG { +\E)[f(E-(, L )r L (E) 

+f(E-^ R )T R (E)]G ( D \E) (35) 

For a full discussion of the actual implementation of these 
expressions see Ref. [lj§|. We stress here we are taking 
into account the electron-electron interactions in the de- 
vice region by an effective mean-field description at the 
level of DFT. Thus the device Hamiltonian is an effective 
Kohn-Sham one-body Hamiltonian. In this approxima- 
tion to the true many-body problem Eqs. [7] and [8] remain 
valid even out of equilibrium and at finite temperature. 

B. Supercell approach 

In the supercell approach we calculate the electronic 
structure of the device region and the electrodes with ab 
initio electronic structure programs for periodic systems 
that use localized basis sets such as Crsytal or Siesta. 
We first define a one-dimensional periodic system con- 
sisting of the device region as the unit cell, as shown in 
Fig. [Tfa). It is crucial that the device part D contains 
a sufficiently large portion of the nanowire electrodes so 
as to guarantee that the electronic structure of the de- 
vice region and thus the Hamiltonian Hd is the same as 
the electronic structure of the device between two semi- 
infinite nanowires. In other words, we seek to connect 
the two leads L and R far enough away from the scatter- 
ing region where the electronic structure has relaxed to 
that of a bulk (i.e., infinite) nanowire. 

In a similar way, the unit cell Hamiltonian matrix 

L /R 

H and hopping matrices between unit cells for left 
and right leads H^ R arc extracted from periodic cal- 
culations of infinite nanowires of finite width [see Fig. 
[7{b,c)]. The lead self-energies £l, £r which describe 
the coupling of the device region D to the semi-infinite 
nanowires L and R in the situation depicted in Fig. [7[d) 
can now be calculated by the Dyson equations (l2~7| and 
(US]). 
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Sequence I 



Sequence II 



Sequence III 



FIG. 7: Illustration of the supercell approach to calculate 
the electronic structure of the device and of the leads: (a) 
One-dimensional periodic system to calculate the electronic 
structure of the device region. (b,c) Infinite nanowires to 
calculate the electronic structure of the left (L) and right (R) 
semi-infinite leads, (d) Sketch of the setup of the physical 
system: The device region (D) is suspended between two semi- 
infinite leads L and R. 



Since the electronic structure of the electrodes has been 
calculated for perfect nanowires, the effect of an eventual 
charge transfer within the device region has not been 
taken into account yet. As pointed out before, the net 
effect of the charge transfer, mostly within the device 
region, on the electronic structure of the bulk electrode 
material outside the device can be taken into account 
by simply shifting the electrochemical potentials (and of 
course the band structure) of the two materials to a com- 
mon electrochemical potential. Hence, the Kohn-Sham 
Hamiltonian of the entire system is also given by Eq. 
(|29|) where now the common electrochemical potential 
jj, is the one obtained from the supercell calculation of 
the device region, and /i° and ^ are the electrochemical 
potentials obtained from the electronic structure calcu- 
lations of the infinite nanowires. The Green's function 
of the device region is given by eq. (|32|) where the self- 
energies Si and S r are now obtained from the Dyson 
equations (|27| and (|28]l but with the energies shifted by 



the electrostatic shifts Al = (J. — /i^ and Ar = f-i 



H! 



( 2 + a l )s; 



(36) 



({z + Al) So - Ho - S^z))" 1 ((Hi - (z + A L ) Si)) 



E r (z) = (Hr - (z + A R ) Si) 

x ((* + A R ) So - Ho - S^z))" 1 (h* 



(37) 
(z + A R )S t 1 



By this procedure we have connected the device region D 
with two semi-infinite nanowires that have the electronic 
structure of bulk, i.e. infinite, nanowires far from the 
device. 




FIG. 8: Sketch of the pyramidal nanocontact geometries used 
in the ab initio quantum transport calculations with Bethe 
lattice electrodes. Each column represents a sequence where 
the amount of bulk electrode material in the device region is 
increased in a disctinct way. 



V. COMPARISON OF ELECTRODE MODELS 

We now compare results obtained with the two elec- 
trode models, i.e., the Bethe lattice and the nanowire 
electrodes and accompanying different implementations 
of the self-consistent procedure. We have chosen for this 
study an archetypal metallic nanocontact model formed 
by two pyramidal tips joined by the apex atoms. For the 
ab initio calculation of the device region and the nanowire 
electrodes we use LDA and the minimal valence basis set 
by Christian and coworkers [8l|. Should one be inter- 
ested in a more quantitative study of the conductance 



of these systems it would be convenient to increase the 
size of the basis set, but this is not the main aim of this 
work. For the Bethe lattice we take the tig ht-binding 
parametrization by Papaconstantopoulos [78| . obtained 
by fitting tight-binding parametrizations to DFT calcula- 
tions. Differences between the different parametrizations 
discussed above are essentially irrelevant. 

In Fig. [8] we show three different sequences of in- 
creasing size for the embedded cluster calculations with 
Bethe lattice electrodes. The narrowest section or con- 
tact atomic structure is maintained throughout the se- 
quence. In sequence I, the pyramidal form is maintained 
for the entire device while it is increased in size. In this 
case the Bethe lattices are only connected to the base 
layers of the pyramids. In sequence II and III, only the 
tips maintain the pyramidal form. The rest of the device 
region arc finite sections of 001 surfaces of the fee crys- 
tal lattice. In each step of a sequence an atomic layer is 
added. The difference between sequences II and III is the 
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18 




FIG. 9: Sketch of the nanocontact geometries and the corre- 
sponding nanowire electrodes used in the ab initio quantum 
transport calculations with nanowire electrodes. 



width of the finite sections of bulk electrode included in 
the device region. 

In Fig. |9]the sequence of model geometries for the case 
of nanowire electrodes is shown. The device region is also 
composed of two pyramidal tips and also contains the 
unit cells used for the computation of the semi infinite 
nanowire self-energies. In each step the section of the 
nanowire is increased. 



Fig. [TOld) shows the transmission functions for Cu 
nanocontacts with the same basic contact geometry as 
before but now with Cu nanowires of finite section serv- 
ing as bulk electrodes instead of the Bethe lattices. Now 
the transmission functions are much spikier than before, 
especially at higher energies. This is somewhat reduced 
by increasing the width of the nanowire electrodes since 
the density of peaks increases and begin to merge. The 
origin of this fine structure in the transmission func- 
tion lies in the well-defined conducting channels in the 
electrodes [Hf. Interestingly the plateau of transmission 
one near the Fermi level is clearly visible. Also the overall 
transmission curves are roughly similar to the ones ob- 
tained before with the Bethe lattice electrodes, at least 
for the bigger nanowires. The computational effort is, 
however, greatly increased in the latter case. 

The relative stability of the T(E) f» 1 plateau near 
the Fermi level with respect to changes in the size and 
specific form of the device region or with respect to the 
electrode model is of course due to the low sensitivity to 
elastic scattering of the very delocalized s-type conduc- 
tion electrons. The variations of the transmission curves 
for the different electrode models and device regions can 
still be attributed to the different interference patterns 
of the conduction electrons. These results are consistent 
with experimental evidence which shows a sharp but nev- 
ertheless finite- width peak in the conductance histogram 
of Cu nanocontacts at 1 Gq. 



B. An sp-type conductor: Al 



A. An s-type conductor: Cu 

First we have studied the relatively simple situation 
of an s-type material, i.e., only s-type electrons are con- 
tributing to the DOS near the Fermi level and hence to 
the zero-bias conductance. Here we consider Cu which is 
a low-resistivity metal frequently used in nano-electronics 
and STM experiments. Figs. |T07 a)-(c) show the trans- 
mission functions (near the Fermi level) calculated with 
Bethe lattice electrodes. The nanocontacts share the 
same basic contact geometry but have different amounts 
of bulk electrode material included in the device region 
according to the different geometry sequences explained 
above and illustrated in Fig. [8] As can be seen the in- 
dividual transmission functions vary for different geome- 
tries within each sequence and also between sequences. 
However, the overall shapes of the individual transmis- 
sion functions are very similar. Near the Fermi level all 
transmission functions feature a plateau around one im- 
plying a zero-bias conductance of approximately 1 Go- 
The origin of this "quantized" conductance lies in the 
single open channel composed of Cu 4s-orbitals which 
is almost perfectly conducting. The number and orbital 
composition of the conducting channels can be done as 
explained in Ref. l82l . 



We now turn to a slightly more complicated situa- 
tion of nanocontacts made from Al which is an sp-type 
conductor. In this case we expect that changes in the 
geometry of the nanoconctact should have a bigger ef- 
fect on the transmission than in the case of a simple s- 
type conductor since p-orbitals contributing now to the 
conductance are more susceptible to elastic scattering 
than s-orbitals due to their directionality. In fact, this 
can be appreciated in the different sequences (see Fig. 
Hip . Now. as the device region increases, the transmis- 
sion curves present larger variability. Furthermore, the 
conductance at the Fermi level changes noticeably, ap- 
proaching a definite value only for large systems where we 
consistently obtain zero-bias conductances below 0.5 Go- 
This result does not seem to be in agreement with ex- 
perimental evidence showing, typically, zero-bias conduc- 
tances between 0.5 and 1 Go for the last conductance 
plateau before breaking [83|. We remind the reader, how- 
ever, that a much more complete statistical analysis is 
needed to draw conclusions in this regard. This anal- 
ysis was carried out in the past [13] and a good agree- 
ment was found between theory and the experimental 
results. A common feature to most transmission curves 
is a pronounced increase above the Fermi level. This 
shoulder or peak originates from the doubly-degenerate 

iPy-channcl while the transmission through the s- and 
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FIG. 10: Transmission functions for Cu nanocontacts with the same tip geometry but for different electrode models, (a)-(c) 
Transmission functions calculated with Bethe lattice electrodes for different amounts of bulk electrode material included into the 
device region according according to the three sequences of geometries shown in Fig. [5] (d) Transmission functions calculated 
with nanowire electrodes of different diameters according to Fig. [9] (the individual transmission curves have been offset by 1 in 
order to distinguish them from each other). 



p z -channels is suppressed [851] . Again the transmission 
functions obtained with the finite-section nanowire elec- 
trodes present more structure [see Fig. Illf d)]. eventually 
converging to the overall shape obtained with the Bethe 
lattice electrodes for large enough section electrodes. 



C. An sd-type conductor: Ni 

As a last example we consider the case of a nanocon- 
tact made out of Ni, an sd material. This case is the more 
complex from the point of view of the electronic structure 
since, in principle, 6 orbitals are expected to contribute 
to the conductance at the Fermi level. For simplicity's 
sake, we consider paramagnetic Ni, where the two spin 
channels contribute equally to the current. A realistic 
theoretical treatment aiming at understanding the avail- 
able experimental results of truely magnetic Ni has been 
presented in the past[H, H?} and will not be repeated 



here. Contrary to what one might expect, the transmis- 
sion at the Fermi level does not depend too much on the 
chosen size of the central region or device when using 
Bethe lattice electrodes, keeping a fairly constant value 
around 3 even for the smallest systems. It is, again, in 
the case of using nanowire electrodes that becomes ap- 
parent that one needs to increase the section of the wires 
before converging to a similar number, at the concommi- 
tant computational cost. The orbital analysis indicates 
that the ,s-channel and two d-channels arc mainly the 
ones responsible for this number. Away from the Fermi 
level the variations among curves are similar to the ones 
in the previous cases, hardly exceeding in magnitude the 
transmission unity. 
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FIG. 11: Transmission functions for Al nanocontacts with the same tip geometry but for different electrode models, (a)-(c) 
Transmission functions calculated with Bethe lattice electrodes for different amounts of bulk electrode material included into 
the device region according to the three sequences of geometries shown in Fig. [S] (d) Transmission functions calculated with 
nanowire electrodes of different diameters according to Fig. [9] (the individual transmission curves have been offset by 1 in order 
to distinguish them from each other). 



VI. DISCUSSION AND CONCLUSIONS 

We have presented a detailed account of the theoret- 
ical and computational treatment of quantum transport 
in nanostructures. While similar analysis have been re- 
ported in the past, ours mainly focuses on implementa- 
tion details usually skipped in the literature, but crucial 
for those interested in developing codes as the two pre- 
sented here. In addition to the well-known pitfalls in the 
use of DFT in transport problems, the way this imple- 
mentation is carried out determines, to a good extent, 
the results or the difficulty in obtaining reliable results 
that can be compared with experiments. We have thus 
made a critical comparison between to archetypal types 
of electrodes, which is a source of discrepancy and con- 
troversy: parametrized vs. ab initio. Without pretend- 
ing our two codes to be representative of all the other 
developed by many groups, we can conclude that the 
use of parametrized electrodes presents two advantages 



with respect to a more faithful description of the elec- 
tronic structure of the electrodes. First, the variability 
between transmission curves is greatly reduced even for 
small devices or central regions when compared to the use 
of nanowire electrodes. Second, the computational cost 
in the calculation of the self-energy in the former case 
can be orders of magnitude smaller than in the latter, 
particularly for large section wires. The use of semiinfi- 
nite wires as electrodes is, nevertheless, essential to prop- 
erly understand scattering in a variety of systems such as 
true semiconducting nanowires, atomic chains [88j. car- 
bon nanotubes, or graphene nanoribbons [8fjj . 
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Appendix A: Representation of operators in 
non-orthogonal basis sets 

The natural definition for the matrix A of a one-body 
operator A in a non-orthogonal basis set (NOBS) {|a)} 
is simply by its matrix elements: 



A = (A a p) = ((a\A\p) 



(Al) 



However, the representation of an operator in a NOBS is 
not that simple: 



i = ^|a)(S- 1 AS- 1 ) Q(3 (/3| 

a,/3 



(A2) 



where S = {S a p) — ((a | /?)) is the overlap matrix. It is 
easy to see that this definition leads results in the matrix 
elements A a p defined above. Then the identity operator 



in the NOBS representation is given by 
I = ^|a>(S- 1 )^</3|, 

which is also easy to proof. 

Now we define a second matrix 

A := S _1 AS _1 , 



(A3) 



(A4) 



which is the matrix that appears above in the represen- 
tation of the operator in a NOBS. 

One should take care when using the representation 
of an operator in a NOBS. For example, the matrix ele- 
ment of an operator between two non-orthogonal orbitals 
a), LS) can be zero, (aL4|/?) = 0, but the correspond- 
ing matrix element of the matrix A does not necessarily 
vanish due to the multiplication with the inverse of the 
overlap matrix on both sides. Thus there is actually a 
non-zero contribution of the two orbitals to the operator 
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although the corresponding matrix element of the oper- 
ator is zero 

Orthogonalizing the basis set by the Lowdin orthogo- 
nalization scheme [9(|, the matrices A and A are trans- 
formed to the matrix A x = f j/j defined in the new 
orthogonal basis set according to: 

A ± = s -i/a AS -i/a = 8+1/2^8+1/2^ (A5) 

Though there are also other orthogonalization schemes, 
the Lowdin scheme is particularly useful in the context of 
quantum chemistry methods based on atomic orbitals as 
the center of the orthogonalizcd orbital remains centered 



on the same atom as the original non-orthogonal orbital. 



Appendix B: Partitioning method 

As explained in Sec. |n]we model the transport prob- 
lem by dividing the system in three parts. Two semi- 
infinite leads (L) and (R) with bulk electronic structure 
are connected to a finite region called device (D). In a 
local basis set the Hamiltonian and the overlap matrix 
of the system are given by ([T]) and (|2|). Dividing the F 
matrix into sub-matrices in a similar manner we obtain 
the following matrix equation: 



z S L - H L 
z Sdl — Hdl 
Orl 



Sld 


— Hld 


zS d 


H D 


Srd 


— Hrd 



Orl 



z S D r - 

^Hr 




G l (s) G ld ( 2 ) G lr (2 
Gdl(z) Go(z) G D r(*, 
Gr L (z) G rd (z) G r (z) 



1l Old L r 
Odl Id Odr 
Orl Ord 1r 



This yields 9 equations for the 9 sub-matrices of the GF 
G. We can resolve this matrix equation columnwise. 
Multiplying all rows of ES — H with the first column 
of G yields three equations for Gl , Gdl and Grl which 
yield: 

G L (z) = (zS L - H L - So+RMr 1 
Gul(z) = g D + R(^)(H DL - 2 S DL )G L (z) 
G RL (z) = g R (z) (Hr D - zSrd) G dl {z) 

Similarly we obtain from multiplication with the second 
column: 

G D (z) = (zSd-Hd-Sl(^)-Sr(^))" 1 
Gld(^) = Sl(z) (H ld - zS LD ) G D (z) 
G RD (z) = Sr(z) (Hrd - 2;S RD ) G d (z) 

And finally from multiplication with the third column, 
we obtain: 

G R (z) = (zSr-Hr-Sd+lW)- 1 
G DR (z) = Ed+l(z) (H D r - zSbr) Gr(z) 
Glr(^) = Eh(z) (H L d - ^Sld) G DR (z) 

We have introduced the Green's functions of the iso- 
lated left and right lead gL and g R and the corresponding 
self-energies El and £r: 

gI » = (zSl-Hl)- 1 

S L (z) ee (H DL -zS DL )g L (z)(H LD -zS LD ) 

gR (z) EE (zSr-Hr)- 1 

Er(z) ee (H D r - zS DR ) g R (z) (Hrd - zS RD ) 



Furthermore, we have defined the Green's function of 
the device plus the left lead only, gD+L, of the device 
plus the right lead only, gD+R, and the corresponding 
self-energies Sd+l and Sd+r each one representing the 
coupling of one of the leads to the device and the other 
lead: 

g D+L (z) ee (zSd-Hd-El(z))- 1 
g D+R (z) ee (zSd-Hd-SrM)- 1 
£ D +r(z) = (Hrd - zSRX>)gD+L(2)(H D R - zS D r) 
S d +l(z) = (Hld - 2;S L D)gD+R(^)(H D L - ^S D l) 

Appendix C: Bethe lattices 

In this appendix we discuss how self-energies for Bethe 
lattices (BL) used to describe the leads are calculated. A 
BL is generated by connecting a site with N nearest- 
neighbors in directions that could be those of a partic- 
ular crystalline lattice. The new N sites are each one 
connected to iV — 1 different sites and so on and so forth. 
The generated lattice has the actual local topology (num- 
ber of neighbors and crystal directions) but has no rings, 
and thus does not describe the long range order charac- 
teristic of real crystals. Let n be a generic site connected 
to one preceding neighbor n — 1 and N — 1 neighbors 
of the following shell (n + i with i = 1,..,N — 1). For 
simplicity's sake, we carry out the derivation for an or- 
thogonal basis. Following App. [Dj the generalization to 
the case of a non-orthogonal basis is straigthforward. 
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Dyson's equation for an arbitrary non-diagonal Green's 
function is 



(EI - H )G„^ 



H, 



E 

= 1 iV- 



(CI) 

where k is an arbitrary site, E the energy, and H 2 J is a 
matrix that incorporates interactions between orbitals at 
sites i and j (bold capital characters are used to denote 
matrices). Ho is a diagonal matrix containing the orbital 
levels and I is the identity matrix. Then, we define a 
transfer matrix as 

Ti-i^Gi-i j = Gjj (C2) 
Multiplying Eq. ()Clj) by the inverse of G n _i. n we obtain, 



(M-Ho)T„ 



E 

=l,...,iV-l 



11- . I* • 



T 



l.n 



(C3) 

Due to the absence of rings the above equation is valid 
for any set of lattice sites, and, thus, solving the BL is 
reduced to a calculation of a few transfer matrices. Note 
that a transfer matrix such as that of Eq. (|C2[) could 
also be defined in a crystalline lattice but, in that case it 
would be useless. 

Eq. (|C3p can be solved iteratively, 



matrix in these directions. To make connection with the 
notation used above note that the vector that joins site 
n — 1 to site n, namely, r n — r n _i would necessarily be 
one of the lattice directions of the set tj. The self-energies 
associated to each direction have to be obtained from the 
following set of 2N coupled self-consistent equations, 



£ Tj = H n [EI H (Sf - S^)]" 1 Hj. 



Sf^ — [EI — Hq — {^i 



SO] _1 Ht 



(C7) 



(C8) 



where i = 1,...,JV and t, = — r,. H Ti is the interatomic 
interaction in the r, direction, and Ex and are the 
sums of the self-energy matrices entering through all the 
Caylcy tree branches attached to an atom and their in- 
verses, respectively, i.e. 



N 



N 



S Ti and 



E^ 



(C9) 



This set of 2N matrix equations has to be solved iter- 
atively. It is straightforward to check that, in cases of 
full symmetry, it reduces to the single equation. The lo- 
cal density of states can be obtained from the diagonal 
Green's function matrix, 



EI — Hn 



E 



=1 N-l 



H; 



(C4) 

If the orbital basis set and the lattice have full symme- 
try (including inversion symmetry) the different transfer 
matrices can be obtained from just a single one through 
appropriate rotations. However this is not always the 
case (see below). 

Before proceeding any further we define self-energies 
that can be (and commonly are) used in place of transfer 
matrices, 



(C5) 



Eq. (|C4[) is then rewritten as, 



H 



n — 1 , n 



EI - H 



E 

i=l N-l 



o 



H 



(C6) 



ei n r 



E ^ 



,N 



(CIO) 



Appendix D: Self-energy of a one-dimensional lead 

Here we will derive the Dyson equation (f2"5| for the 
calculation of the self-energy of the semi-infinite right 
lead. The derivation of the Dyson equation for the left 
lead l|27p goes in a completely analogous way. 

The Hamiltonian matrix Hr of the (isolated) semi- 
infinite right electrode is defined in eq. (j2T))) as: 



Hr 



/H Hi 
H\ H Hi 



Hi H Hi 



V o 



\ 



/ 



(Dl) 



where we have made use of the general property and the overlap matrix is given in eq. (22) as: 



H 



H 



n,n-l — *-*-n-l,n- 

As discussed hereafter, in a general case of no symme- 
try this would be a set of N coupled equations (2iV if 
there is no inversion symmetry). Symmetry can be bro- 
ken due to either the spatial atomic arrangement, the or- 
bitals on the atoms that occupy each lattice site, or both. 
When no symmetry exists, the following procedure has 
to be followed to obtain the self-energy in an arbitrary 
direction. The method is valid for any basis set or lattice. 
Let 71 be the N nearest-neighbor directions of the lattice 
we are interested in and H T . the interatomic interaction 



/So Si ^ 

Si So Si 
Si Sq Si 



V o 



(D2) 



/ 



To obtain the self-energy of the lead we have to calculate 
the GF of the lead from its defining equation: 



(zS R - H R ) gR (z) = 1 



(D3) 
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In the same way as the Hamiltonian and the overlap matrix we subdivide the GF matrix gR into sub-matrices 
corresponding to the unit cells of the lead. Now the above equation for the right lead's GF reads: 



/ zS - H zSi - Hi 
zSt - U] zS - H zSi - Hi 



gl,l gl,2 
g2,l g2,2 



/I o 
1 



(D4) 



As explained in Sec lIII Bl it suffices to calculate the "surface" GF, i.e. gi^. From multiplication of the 1st, the 2nd 
and so on until the n-th line of (zSr — Hr) with the 1st column of gR(z) we get the following chain of equations: 



(zS o -Ho)gi,i(20 + («Si-Hi)g2,i(«) = 1 
(zS t 1 -H t 1 )g 1 , 1 (z) + (zSo-Ho)g 2 a(z) + (zS 1 -H 1 )g 3 , 1 (z) = 



(D5) 
(D6) 



(zS t 1 -H t 1 )g„_ 1 , 1 (z) + (zS -H )g„,i(z) + (zSi-H 1 )g n+1 , 1 (z) = 
For n > 1 the equations for determining g„.i(z) all have the same structure: 

(zSo-H )g„,i(z) = (H t 1 -zS t 1 )g n _ lil (z) + (H 1 -zS 1 )g n+1 , 1 (z). 



(D7) 



(D8) 



We define a transfer matrix for n > 1 by: 

T„_i,„(«)g n _i,i(z) = g n ,i(z). (D9) 

The transfer matrix thus transfers information from site 
n — 1 to site n of the lead, i.e. from the left to the right. 
Multiplying Eq. (|D8[) by (gn-i,i) we obtain: 

(zSo-H )T„_i, n (z) = (D10) 
(Hj - zS\) + (Hi - zSi) T„,„ +1 (z) T n _ 1>n (z) 

Reordering we obtain the following iterative equation for 
the transfer matrices: 

T n _i,„(z) = (Dll) 
(zS H (Hi - zSr) T^+^z))- 1 (Hj - zS\) 

Since the electrode is semi-infinite it looks the same from 
each unit cell when looking to the right. Thus a given 
g n — results always in the same g n i independent of 
n. Thus the transfer matrix must be independent of n: 
T„_x n {z) = T(z), and Eq. (|D11|) allows to determine 
the T(z) sclf-consistcntly. 

We define the self-energy as S(z) := (Hi — zS\) T(z), 
and obtain the Dyson equation for the self-energy: 

S(z) = (Hi - zSi) (zS Ho S(z))- 1 (H+ - zS\) 

(D12) 

We will now see that this self-energy is indeed identi- 
cal to the one defined for the right lead in eq. (|28|). 



i.e. S(z) = H r (E). By plugging in the definition of the 
transfer matrix, eq. (|D9[) . into eq. (|D5[) for determining 
the surface GF, gi^ we find: 

(zSo-Ho)g 1 ,i(^) + (^Si-Hi)T(z)gia(z) = 1 
Plugging in the definition of the self-energy we obtain: 

=> (zS -Ho + S(z))g M (z) = 1 (D13) 

Thus we obtain for the surface GF of the right lead: 

gi.!(z) = (zS -H + S(z))- 1 . (D14) 

And vice- versa the self-energy can be expressed in terms 
of the surface GF: 

S(z) = (Hi - zSi) gi,i(z) (Hi - zS\). (D15) 

This proofs that the self-energy S(z) defined above in 
terms of the transfer matrix is identical to the self-energy 
S r (z) defined earlier in Sec lIII Bl so that the self-energy 
S r (z) can be calculated iteratively by the Dyson equa- 
tion pgj) . 

The proof for the left lead runs completely analogously. 
The surface GF of the left lead is now: 

g-i,-i(z) = (zSo-Ho + E^z))- 1 . (D16) 
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